
#include "glodef.h"
#include "errlog.h"
#include "vecs.h"
#include "grid3d.h"
#include "g3dutils.h"
#include "f3io.h"

#include "soah_common.h"
#include "glovar.h"

#include "ns_csplit.h"

#include <stdio.h>
#include <string.h>

double current_adj_solution_norm[];
double time_inf_adj_solution_norm[];
char adj_nrms_fname[FILENAME_MAX] = "adj_norms.dat";


static void make_adv_u(void);
static void adjoint_norms(void);
static void print_inf_adj_norms(void);

void solve_periodic_adjoint_fvp(const vector *v_finl, scalar *q_finl)
{
    
    curr_ti = total_ti;
    move_grid();
    load_state_time(state_dir());
    
    PrintLogLine(78, '=');
    PrintLog("Time %-10f\n", (double) TIME);
    
    /* make final adjoint pressure compatible with v_final and u */
#if defined(CS_INITIAL_PRESSURE)
    {
        FinalAdjointPressure_cs(solver_data, v_finl, u, NULL, NULL, NULL, q_finl);
        make_scal_mean_0(grid, elemP, q_finl);
    }
#endif
    
    copy_vec(v_finl, v);
    copy_scal(q_finl, q);
    clean_scal(qhi);
    
    if(!quiet)
        printf("\rT = %6.3f; ", (double)(TIME));
    adjoint_norms();
    if(!quiet)
        fflush(stdout);
    
    create_dir(adjoint_dir());
    save_adjoint_time(adjoint_dir());
    
    for(; curr_ti > 0;)
    {
        curr_ti--;
        
        clock_t tm = clock();
        
        PrintLogLine(78, '=');
        PrintLog("Time %-10f\n", (double) TIME);
        
        load_state_time(state_dir());
        
        if(!flags.moving_grid)
        {
            copy_vec(u, adv_u);
        }
        else
        {
            move_grid();
            make_adv_u();
        }
        copy_vec(v, v0);
        copy_scal(q, q0);
        copy_scal(qhi, qhi0);
        
        /* do a timestep */
        TimeStep_adjoint_cs(solver_data, v, qhi, q, adv_u, u, NULL, NULL, v0, qhi0, q0, FALSE);
        /* the solver may set p to 0 in one point.  To compute l2 norms we
         * need p to have mean 0, so
         */
        make_scal_mean_0(grid, elemP, q);
        
        save_adjoint_time(adjoint_dir());
        
        if(!quiet)
            printf("\rT = %6.3f; ", (double)(TIME));
        adjoint_norms();
        if(!quiet)
            fflush(stdout);
        
        PrintLog("CPU time for this timestep: %-6.3fsec.\n",
                (clock() - tm)/(double) CLOCKS_PER_SEC);
    }
    
}


/****************************************************************************/
/****************************************************************************/
/****************************************************************************/

void solve_adjoint_fvp(const vector *v_finl, scalar *q_finl,
        vector *vtT_bc)
{
    unsigned i;
    double adj_norms[4];
    static vector _src, _vbc;
    vector * const src = &_src;
    vector * const vbc = &_vbc;
    
    (void) vtT_bc;
    
    alloc_vec(src, v->N);
    alloc_vec(vbc, v->N);    
    
    curr_ti = total_ti;
    move_grid();
    
    if(flags.conv_adjoint)
    {
        interpolate_vector_lagr(grid, elemU, TIME, exct_vel, u);
        interpolate_scalar_lagr(grid, elemP, TIME, exct_prs, p);
        make_scal_mean_0(grid, elemP, p);        
    }
    else if(dump_bin > 0)
        load_state_time(state_dir());
    
    PrintLogLine(78, '=');
    PrintLog("Time %-10f\n", (double) TIME);
    
    /* make final adjoint pressure compatible with v_final and u */
#if defined(CS_INITIAL_PRESSURE)
    {
        interpolate_vector_lagr(grid, elemU, TIME, srce_trm, src);
        for(i=0;i<3*src->N;src->X[i++]*=-1);
        FinalAdjointPressure_cs(solver_data, v_finl, u, vtT_bc, NULL, NULL, q_finl);
        make_scal_mean_0(grid, elemP, q_finl);
    }
#endif
    
   
    copy_vec(v_finl, v);
    copy_scal(q_finl, q);
    clean_scal(qhi);

    if(!quiet)
        printf("\rT = %6.3f; ", (double)(TIME));

    memset(time_inf_adj_solution_norm, 0, sizeof(time_inf_adj_solution_norm));
    adjoint_norms();

    if(dump_bin > 0)
    {
        create_dir(adjoint_dir());
        save_adjoint_time(adjoint_dir());
    }
    for(; curr_ti > 0;)
    {
        curr_ti--;
        
        clock_t tm = clock();
        
        PrintLogLine(78, '=');
        PrintLog("Time %-10f\n", (double) TIME);
        
        if(flags.conv_adjoint)
        {
            interpolate_vector_lagr(grid, elemU, TIME, exct_vel, u);
            interpolate_scalar_lagr(grid, elemP, TIME, exct_prs, p);
            make_scal_mean_0(grid, elemP, p);        
        }
        else if(dump_bin > 0)
            load_state_time(state_dir());
        
        if(!flags.moving_grid)
        {
            copy_vec(u, adv_u);
        }
        else
        {
            move_grid();
            make_adv_u();
        }
        copy_vec(v, v0);
        copy_scal(q, q0);
        copy_scal(qhi, qhi0);
        
        test_vector_field(grid, elemU, TIME, adj_srce_trm, src);
        apply_bc_dir_vector(grid, elemU, TIME, adj_bdry_vel, vbc);
        
        /* make a timestep */
        TimeStep_adjoint_cs(solver_data, v, qhi, q, adv_u, u, src, vbc, v0, qhi0, q0, FALSE);
        make_scal_mean_0(grid, elemP, q);
        
        if((dump_bin > 0) && (curr_ti % dump_bin)==0)
            save_adjoint_time(adjoint_dir());
        
        if(!quiet)
            printf("\rT = %6.3f; ", (double)(TIME));
        
        adjoint_norms();

        if(!quiet)
            fflush(stdout);
        
        PrintLog("CPU time for this timestep: %-6.3fsec.\n",
                (clock() - tm)/(double) CLOCKS_PER_SEC);
    }
    print_inf_adj_norms();
    free_vec(src);
    free_vec(vbc);    
}



/****************************************************************************/
/****************************************************************************/
/****************************************************************************/
/****************************************************************************/

void solve_adjoint_11_fvp(const vector *v_finl, scalar *q_finl,
        vector *vtT_bc)
{
    unsigned i;
    double adj_norms[4];
    vector *vbc = u;
    static vector _src;
    vector *src = &_src;
    
    (void) vtT_bc;
    
    alloc_vec(src, u->N);
    
    curr_ti = total_ti;
    move_grid();
    load_state_time(state_dir());
    
    PrintLogLine(78, '=');
    PrintLog("Time %-10f\n", (double) TIME);
    
    /* make final adjoint pressure compatible with v_final and u */
#if defined(CS_INITIAL_PRESSURE)
    {
        interpolate_vector_lagr(grid, elemU, TIME, srce_trm, src);
        for(i=0;i<3*src->N;src->X[i++]*=-1);
        FinalAdjointPressure_cs(solver_data, v_finl, u, vtT_bc, NULL, NULL, q_finl);
        make_scal_mean_0(grid, elemP, q_finl);
    }
#endif
    
    clean_scal(q_finl);
    
    copy_vec(v_finl, v);
    copy_scal(q_finl, q);
    clean_scal(qhi);
    
    interpolate_vector_lagr(grid, elemU, TIME, exct_vel, eu);
    interpolate_scalar_lagr(grid, elemP, TIME, exct_prs, ep);
    make_scal_mean_0(grid, elemP, ep);
    adj_norms[NORM_U_L2] = norm_vec_diff_l2(grid, elemU, eu, v);
    adj_norms[NORM_U_H1] = norm_vec_diff_h1(grid, elemU, eu, v);
    adj_norms[NORM_P_L2] = norm_scal_diff_l2(grid, elemP, ep, q);
    adj_norms[NORM_P_H1] = norm_scal_diff_h1(grid, elemP, ep, q);
    
    PrintLogNorm("v-eu", "L2", adj_norms[NORM_U_L2]);
    PrintLogNorm("q-ep", "L2", adj_norms[NORM_P_L2]);
    PrintLogNorm("v-eu", "H1", adj_norms[NORM_U_H1]);
    PrintLogNorm("q-ep", "H1", adj_norms[NORM_P_H1]);
    
    create_dir(adjoint_dir());
    save_adjoint_time(adjoint_dir());
    
    for(; curr_ti > 0;)
    {
        curr_ti--;
        
        clock_t tm = clock();
        
        PrintLogLine(78, '=');
        PrintLog("Time %-10f\n", (double) TIME);
        
        load_state_time(state_dir());
        
        if(!flags.moving_grid)
        {
            copy_vec(u, adv_u);
        }
        else
        {
            move_grid();
            make_adv_u();
        }
        copy_vec(v, v0);
        copy_scal(q, q0);
        copy_scal(qhi, qhi0);
        
        
        /*
         * vector *vbc = u;
         * vector *src = NULL;
         */
        
        interpolate_vector_lagr(grid, elemU, TIME, srce_trm, src);
        for(i=0;i<3*src->N;src->X[i++]*=-1);
        
        /* do a timestep */
        TimeStep_adjoint_cs(solver_data, v, qhi, q, adv_u, u, NULL, vbc, v0, qhi0, q0, TRUE);
        /* the solver may set p to 0 in one point.  To compute l2 norms we
         * need p to have mean 0, so
         */
        make_scal_mean_0(grid, elemP, q);
        
        save_adjoint_time(adjoint_dir());
        
        if(!quiet)
            printf("\rT = %6.3f; ", (double)(TIME));
        
        interpolate_vector_lagr(grid, elemU, TIME, exct_vel, eu);
        interpolate_scalar_lagr(grid, elemP, TIME, exct_prs, ep);
        make_scal_mean_0(grid, elemP, ep);
        adj_norms[NORM_U_L2] = norm_vec_diff_l2(grid, elemU, eu, v);
        adj_norms[NORM_U_H1] = norm_vec_diff_h1(grid, elemU, eu, v);
        adj_norms[NORM_P_L2] = norm_scal_diff_l2(grid, elemP, ep, q);
        adj_norms[NORM_P_H1] = norm_scal_diff_h1(grid, elemP, ep, q);
        
        PrintLogNorm("v-eu", "L2", adj_norms[NORM_U_L2]);
        PrintLogNorm("q-ep", "L2", adj_norms[NORM_P_L2]);
        PrintLogNorm("v-eu", "H1", adj_norms[NORM_U_H1]);
        PrintLogNorm("q-ep", "H1", adj_norms[NORM_P_H1]);
        if(!quiet)
            printf("||v-eu||_L2 = %10.4e", adj_norms[NORM_U_L2]);
        
        if(!quiet)
            fflush(stdout);
        
        PrintLog("CPU time for this timestep: %-6.3fsec.\n",
                (clock() - tm)/(double) CLOCKS_PER_SEC);
    }
    free_vec(src);
}


/****************************************************************************/
/****************************************************************************/
/****************************************************************************/
/****************************************************************************/
/****************************************************************************/

/*  #include "fe_bubble.h" */
void solve_shape_adjoint_fvp(const vector *v_finl, scalar *q_finl)
{
    unsigned i;
    vector _src, *src = &_src,
            _grad_u, *grad_u = &_grad_u,
            _curl_u, *curl_u = &_curl_u;
    
    alloc_vec(src, u->N);
    alloc_vec(grad_u, u->N);
    alloc_vec(curl_u, u->N);
    
    curr_ti = total_ti;
    move_grid();
    load_state_time(state_dir());
    
    PrintLogLine(78, '=');
    PrintLog("Time %-10f\n", (double) TIME);
    
    /* make final adjoint pressure compatible with v_finl and u */
#if defined(CS_INITIAL_PRESSURE)
    {
        surf_test_n_cross_curl_bub_in_grad_p1(grid, u, q);
        for(i=0;i<q->N;i++) q->F[i] = (alpha+beta) * q->F[i];
        
        FinalAdjointPressure_cs(solver_data, v_finl, u, NULL, NULL, q, q_finl);
        make_scal_mean_0(grid, elemP, q_finl);
    }
#endif
    
    copy_vec(v_finl, v);
    copy_scal(q_finl, q);
    clean_scal(qhi);
    
    if(!quiet)
        printf("\rT = %6.3f; ", (double)(TIME));
    adjoint_norms();
    if(!quiet)
        fflush(stdout);
    
    create_dir(adjoint_dir());
    save_adjoint_time(adjoint_dir());
    
    for(; curr_ti > 0;)
    {
        curr_ti--;
        
        clock_t tm = clock();
        
        PrintLogLine(78, '=');
        PrintLog("Time %-10f\n", (double) TIME);
        
        copy_vec(u, u0);
        copy_scal(p, p0);
        load_state_time(state_dir());
        
        if(!flags.moving_grid)
        {
            copy_vec(u, adv_u);
        }
        else
        {
            move_grid();
            make_adv_u();
        }
        copy_vec(v, v0);
        copy_scal(q, q0);
        copy_scal(qhi, qhi0);
        
        /* compile the right-hand-sode of the adjoint equation */
        test_grad_e1_in_grad_e2_vec(grid, elemU, u, elemU, grad_u);
        test_curl_e1_in_curl_e2_vec(grid, elemU, u, elemU, curl_u);
        for(i=0;i<3*u->N;i++)
            src->X[i] = alpha*grad_u->X[i] + beta*curl_u->X[i];
        
        /* do a timestep */
        TimeStep_adjoint_cs(solver_data, v, qhi, q, adv_u, u, src, NULL, v0, qhi0, q0, FALSE);
        /* the solver may set p to 0 in one point.  To compute l2 norms we
         * need p to have mean 0, so
         */
        make_scal_mean_0(grid, elemP, q);
        
        save_adjoint_time(adjoint_dir());
        
        if(!quiet)
            printf("\rT = %6.3f; ", (double)(TIME));
        adjoint_norms();
        if(!quiet)
            fflush(stdout);
        
        PrintLog("CPU time for this timestep: %-6.3fsec.\n",
                (clock() - tm)/(double) CLOCKS_PER_SEC);
    }
    
    free_vec(src);
    free_vec(grad_u);
    free_vec(curl_u);
}

/****************************************************************************/

static void make_adv_u(void)
{
    unsigned i;
    vector *ndisp;
    
    assert(elemU==fe_bub || elemU==fe_bub_vec);
    ndisp = get_ndisp(alegrid);
    assert(ndisp->N==grid->NN);
    
    copy_vec(u, adv_u);
    for(i=0;i<grid->NN;i++)
    {
        adv_u->X[i] += ndisp->X[i] / TStep;
        adv_u->Y[i] += ndisp->Y[i] / TStep;
        adv_u->Z[i] += ndisp->Z[i] / TStep;
    }
}

/****************************************************************************/

static void adjoint_norms(void)
{
    interpolate_vector_lagr(grid, elemU, TIME, adj_exct_vel, ev);
    interpolate_scalar_lagr(grid, elemP, TIME, adj_exct_prs, eq);
    make_scal_mean_0(grid, elemP, eq);
    
    current_adj_solution_norm[NORM_DUDT] =
            norm_vec_diff_l2(grid, elemU, v, v0) / TStep;
    
    current_adj_solution_norm[NORM_U_L2] = norm_vec_diff_l2(grid, elemU, v, ev);
    current_adj_solution_norm[NORM_U_H1] = norm_vec_diff_h1(grid, elemU, v, ev);
    current_adj_solution_norm[NORM_P_L2] = norm_scal_diff_l2(grid, elemP, q, eq);
    current_adj_solution_norm[NORM_P_H1] = norm_scal_diff_h1(grid, elemP, q, eq);
    
    for(int i=0;i<4;i++)
        time_inf_adj_solution_norm[i] = fmax(time_inf_adj_solution_norm[i],
                current_adj_solution_norm[i]);
    
    if(adj_exct_vel != NULL)
    {
        PrintLogNorm("v-ev", "L2", current_adj_solution_norm[NORM_U_L2]);
        PrintLogNorm("q-eq", "L2", current_adj_solution_norm[NORM_P_L2]);
        PrintLogNorm("v-ev", "H1", current_adj_solution_norm[NORM_U_H1]);
        PrintLogNorm("q-eq", "H1", current_adj_solution_norm[NORM_P_H1]);
        if(!quiet)
            printf("||v-ev||_L2 = %10.4e", current_adj_solution_norm[NORM_U_L2]);
    }
    else
    {
        PrintLogNorm("v", "L2", current_adj_solution_norm[NORM_U_L2]);
        PrintLogNorm("q", "L2", current_adj_solution_norm[NORM_P_L2]);
        PrintLogNorm("v", "H1", current_adj_solution_norm[NORM_U_H1]);
        PrintLogNorm("q", "H1", current_adj_solution_norm[NORM_P_H1]);
        if(flags.stationary)
        {
            PrintLogNorm("dvdt", "L2", current_adj_solution_norm[NORM_DUDT]);
            if(!quiet)
                printf("||dvdt||_L2 = %10.4e", current_adj_solution_norm[NORM_DUDT]);
        }
        else if(!quiet)
            printf("||  v ||_L2 = %10.4e", current_adj_solution_norm[NORM_U_L2]);
        return;
    }
}

static void print_inf_adj_norms(void)
{
    
    PrintLogLine(1, ' ');
    PrintLogLine(78, '=');
    if(adj_exct_vel!=NULL)
    {
        PrintLogNorm("v-ev", "L2Linf", time_inf_adj_solution_norm[NORM_U_L2]);
        PrintLogNorm("q-eq", "L2Linf", time_inf_adj_solution_norm[NORM_P_L2]);
        PrintLogNorm("v-ev", "H1Linf", time_inf_adj_solution_norm[NORM_U_H1]);
        PrintLogNorm("q-eq", "H1Linf", time_inf_adj_solution_norm[NORM_P_H1]);
    }
    else
    {
        PrintLogNorm("  v ", "L2Linf", time_inf_adj_solution_norm[NORM_U_L2]);
        PrintLogNorm("  q ", "L2Linf", time_inf_adj_solution_norm[NORM_P_L2]);
        PrintLogNorm("  v ", "H1Linf", time_inf_adj_solution_norm[NORM_U_H1]);
        PrintLogNorm("  q ", "H1Linf", time_inf_adj_solution_norm[NORM_P_H1]);
    }
    PrintLogLine(78, '=');
    
    if(flags.conv_adjoint)
    {
        FILE *F = fopen(adj_nrms_fname, "a+t");
        if(F == NULL)
            return;
        fseek(F, 0, SEEK_SET);
        if(fgetc(F)==EOF)
        {
            fprintf(F, "%-12s    %-20s    %-14s  %-14s  %-14s  %-14s\n",
                    "dT", "grid", "||v-ev||L2linf", "||q-eq||L2linf",
                    "||v-ev||H1linf", "||q-eq||H1linf");
        }
        else
            fseek(F, 0, SEEK_END);
        /* fprintf(F, "%-12g    %-12g", (double) TStep, (double) GStep); */
        fprintf(F, "%-12g    %-20.20s", (double) TStep, grid->gridname);
        fprintf(F, "    %12e", time_inf_adj_solution_norm[NORM_U_L2]);
        fprintf(F, "    %12e", time_inf_adj_solution_norm[NORM_P_L2]);
        fprintf(F, "    %12e", time_inf_adj_solution_norm[NORM_U_H1]);
        fprintf(F, "    %12e", time_inf_adj_solution_norm[NORM_P_H1]);
        fprintf(F, "\n");
        fclose(F);
    }
}

